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We theoretically analyze a quasi-two-dimensional system of fermionic polar molecules in a harmonic trans- 
verse confining potential. The renormalized energy bands are calculated by solving the Hartree-Fock equa- 
tion numerically for various trap and dipolar interaction strengths. The inter-subband excitations of the sys- 
tem are studied in the conserving time-dependent Hartree-Fock (TDHF) approximation from the perspective 
of lattice modulation spectroscopy experiments. We find that the excitation spectrum consists of both inter- 
subband particle-hole excitation continuums and anti-bound excitons, arising from the anisotropic nature of 
dipolar interactions. The excitonic modes capture the majority of the spectral weight. We also evaluate the 
inter-subband transition rates in order to investigate the nature of the excitonic modes and find that they are anti- 
bound states formed from particle-hole excitations arising from several subbands. Our results indicate that the 
excitonic effects are present for interaction strengths and temperatures accessible in current experiments with 
polar molecules. 



I. INTRODUCTION 

In the past decade, much of the experimental and the- 
oretical progress in the field of ultracold atoms (DO are 
motivated by the prospect of realizing novel strongly corre- 
lated many-body states of matter. In particular, experimental 
realization of quantum degenerate gases of fermionic polar 
molecules has witnessed a very rapid progress. By associa- 
tion of atoms via a Feshbach resonance to form deeply bound 
ultracold molecules El El, a nearly degenerate gas of KRb po- 
lar molecules in their rotational and vibrational ground state 
has been recently realized Il6l fl0l . The molecules can be po- 
larized by applying a d.c. electric field, resulting in strong 
dipole-dipole inter-molecular interactions. 

A strikingly new feature of systems of fermionic polar 
molecules is the anisotropy of dipolar interactions, making 
them unparalleled among the traditional condensed matter 
systems. Thus, experiments with polar molecules go be- 
yond quantum simulation of effective theories motivated by 
electronic systems and aim at exploring a genuinely new do- 
main of many-body quantum behavior, unique to dipolar in- 
teractions. Dipolar interactions can be utilized to generate 
long-range interactions of arbitrary shape using microwave 
fields [11], simulate exotic spin Hamiltonians [FPU FPU and 
are theoretically predicted to give rise to numerous inter- 
esting collective phenomena such as roton softening lfl4l - 
[T6l . supersolidity lfT7l42Tl . p-wave superfluidity 1221 . emer- 
gence of artificial photons [23 ], bilayer quantum phase transi- 
tions l24l . multi-layer self-assembled chains |25| for bosonic 
molecules, dimerization and inter-layer pairing |26, 27 1, spon- 
taneous inter-layer coherence ||28ll , itinerant ferroelectric - 
ity [29 1, anisotropic Fermi liquid theory and anisotropic sound 
modes I3014331 . fractional quantum Hall effect P4l . Wigner 
crystallization l35l . density-wave and striped order l36l [371 . 
biaxial nematic phase |38|, topological superfluidity ||39l and 
Z2 topological phase [40|, just to mention a few. 

Despite the theoretical prediction of numerous exotic quan- 
tum many-body phenomena in polar molecules, realization 
and observation of many of these novel predictions are still 
an experimental challenge. At the time this paper was written, 



the coldest gas of fermionic polar molecules has been realized 
with KRb molecules at a temperature of lATp llrJ UTOl , where 
Tp is the Fermi temperature. The majority of the mentioned 
quantum phenomena require strong suppression of thermal 
fluctuations, i.e. strong degeneracy condition (T <C Tp). 

A major obstacle towards further evaporative cooling of 
a large class of bi-alkali polar molecules (KRb, LiNa, LiK, 
LiRb, and LiCs) is the existence of an energetically allowed 
two-body chemical reaction channel [41], resulting in signifi- 
cant molecule losses in two-body scatterings. In a low temper- 
ature gas composed of a single hyperfine state, Fermi statis- 
tics blocks scatterings in the s-wave channel and the major- 
ity of scatterings take place through the p-wave channel. In 
unstructured three-dimensional clouds, the attractive head-to- 
tail dipolar interactions soften the p-wave centrifugal barrier 
and increase the cross section of reactive collisions. However, 
the rate of chemical reactions can be effectively suppressed 
by loading the gas into a one-dimensional optical lattice (or 
trap) and aligning the dipoles perpendicular to the formed 
pancakes. In such geometries, the incidence of head-to-tail 
scatterings is effectively suppressed due to the transverse con- 
finement of the gas on one hand, and reinforcement of the p- 
wave barrier due to repulsive side-by-side dipolar interactions 
on the other hand J5] 021 |43j . Therefore, the preferred geom- 
etry to study reactive polar molecules is in tightly confined 
two-dimensional layers. 

Recently, it has been shown that the suppression of chemi- 
cal reactions of reactive fermionic molecules in confined ge- 
ometries remains effective even in the quasi-two-dimensional 
limit, i.e. where not only the lowest, but also the first few ex- 
cited subbands (transverse modes) are populated [43|. On the 
other hand, occupation of higher subbands does not impose 
any difficulty on experiments with non-reactive molecules 
such as NaK, NaRb, NaCs, KCs, and RbCs ED. The pos- 
sibility of going beyond the single-subband limit opens a new 
window towards experimental and theoretical exploration of 
many-body physics of quasi-two-dimensional fermions with 
anisotropic interactions. 

Besides being highly anisotropic, fermionic polar 
molecules and electronic systems differ in another im- 
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portant way: in contrast to the Coulomb interactions, 
dipole-dipole interactions dominate the kinetic energy in 
the high density limit. In d dimensions, dipolar interaction 
scales as r?l d (n is the density) while the kinetic energy 
scales as n 2 / d , implying that their ratio scales as n 1 / d , i.e. 
the interactions are more pronounced at higher densities. We 
define the following dimensionless quantity as a measure of 
dipolar interaction in two dimensions: 



mD 2 n 1 / 2 



(1) 



where D is the electric dipole moment of a single molecule. 

The main goal of this paper is studying the inter-subband 
collective modes of quasi-two-dimensional polar molecules 
and to propose experimental signatures of such phenomena. 
At the time this paper was written, the dipolar interaction 
strengths accessible in the experiments belong to the weakly 
interacting regime (e.g. w 0.05 in experiments of the group 
at JILA) and therefore, we restrict our analysis to the same 
regime. For small r^, the normal liquid phase is expected to 
be the stable phase of the system. 

In the first part of this paper, we study the effect of dipo- 
lar interactions on the single-particle energy dispersions by 
solving the Hartree-Fock equation. In order to simplify our 
analysis, we focus on a single pancake and neglect compli- 
cations such as inter-layer tunneling and attractive inter-layer 
interactions in this study. This approximation is relevant to a 
well-separated stack of pancakes as well. The normal phase 
and the collective modes of fermionic polar molecules in the 
strictly two-dimensional case (single-subband limit) has been 
recently studied by different authors lt22l l33l l36l . Therefore, 
we do not discuss the intra-subband collective modes here and 
instead, focus on genuinely quasi-two-dimensional phenom- 
ena. 

In the second part of this paper, we explore the inter- 
subband excitation spectrum from the lattice modulation 
spectroscopy perspective, an experimental technique origi- 
nally developed for studying strongly correlated phases of 
cold atomic systems I RTIFT 81 . We predict the energy absorp- 
tion and inter-subband transition rates for an exponentially 
switched-on a.c. lattice modulation pulse in the conserving 
time-dependent Hartree-Fock (TDHF) approximation ll49l . 
We show that the excitation spectrum consists of both 
inter-subband particle-hole (p-h) excitation continuums and 
anti-bound excitons, arising from the anisotropic nature of 
dipolar interactions. We show that these many-body effects 
persist for weak interactions and temperatures of the order 
of the Fermi temperature, making their observation feasible 
in current experiments. Finally, we study the inter-subband 
particle transition rates in order to shed light on the nature of 
excitons. We find that in contrast to p-h excitation continuums 
which are associated to certain inter-subband transitions, the 
excitons are composite anti-bound states formed from p-h 
excitations arising from several subbands. 

Before we embark on the formal development of the out- 
lined program, we find it worthwhile to summarize the ob- 
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FIG. 1. A comparison between quasi-two-dimensional fermions in 
the absence (left column) and presence (right column) of interac- 
tions, (a) and (b) show a schematic plot of the band structure. The 
blue and red arrows indicate p-h excitations from to 2, and 1 to 3 
subbands respectively, at different regions in the Fermi sea. (c) and 
(d) show a schematic representation of excited states probed in mod- 
ulation spectroscopy experiments. The labels, — > 2 and 1—^3, 
indicate the inter-subband transition associated to the excitation. The 
blue and red intervals denote p-h continuums and the red lines de- 
note excitonic modes, (e) and (f) show the energy absorption rate, E, 
and inter-subband transition rates, i?o^2 and i?i_>3, as a function of 
modulation frequency, o; mo d. The blue and red shaded continuums 
correspond to — > 2 and 1 — > 3 inter-subband transitions respec- 
tively. The red spike denotes an exciton. The kink in the middle of 
the red shaded continuum is indicating of an exciton lying inside that 
continuum (see plots c and d). Refer to the text for more details. 



tained results in a schematic way. Fig.[T]shows a side by side 
comparison between non-interacting and interacting quasi- 
two-dimensional fermions with dipolar interactions. A typical 
plot of the band structure is shown in Fig. [T^ and |TJ?. In the 
absence of interactions, however, the single-particle energy 
dispersions are quadratic and have a constant energy separa- 
tion of hcu tlap , where w tra p is the trap frequency (we neglect 
the anharmonicity of the trap in order to simplify the analy- 
sis). In the presence of interactions, the single-particle energy 
dispersions are renormalized and no longer remain quadratic. 

We will see later that modulation of the optical lattice stim- 
ulates excitation of particles from populated subbands in equi- 
librium to their second next subband (0 — > 2, 1 —> 3, etc). 
The excitation energy of such processes can be directly mea- 
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sured in the experiments. In the absence of interactions, all 
of the elementary inter-subband excitations have a constant 
energy cost of 2hu trap . The interactions modify the degen- 
erate excited states dramatically, transforming the p-h excita- 
tion energies from their constant value of 2hu trap to a collec- 
tion of continuums and discrete collective modes (excitons), 
as shown in Fig.[JJ: and[TJl. Also, the excitonic modes capture 
the majority of the spectral weight, leaving behind a small 
trace of p-h continuums. 

The mechanism of energy absorption in lattice modulation 
experiments can be explored further by studying the inter- 
subband transition rates. If the modulation frequency lies 
within a certain inter-subband p-h excitation continuum, we 
expect to observe transitions only between the two involved 
subbands, leaving the population of other subbands unaf- 
fected. For example, the blue shaded continuum in the energy 
absorption spectrum in Fig. [TJl and[ljis associated to — > 2 
excitations and consequently, the continuum only appears in 
the plot of -Ro->2 and is absent in the plot of Ri^3 vs. mod- 
ulation frequency, where R a ^/3 denotes the net current from 
a'th to /3'fh Hartree-Fock subband. 

If the modulation frequency corresponds to an excitonic 
mode, we expect to observe transition currents between sev- 
eral pairs of subbands, given that the excitons are generally 
mixtures of p-h excitations of different subbands. For exam- 
ple, the red spikes in the energy absorption spectrum, which 
denotes an exciton, is visible in both i?o^2 and i?i_j.3 plots. 

Finally, if an exciton lies inside a continuum, it will be 
damped and yield a broadened peak in the energy absorption 
spectrum. In this case, although the continuum is associated 
to a certain inter-subband transition per se, we expect to 
observe its trace in other inter-subband transition rate plots 
due to mixing with the excitonic mode. This effect is 
schematically shown in plots Fig. [TJl and[T]\ The red shaded 
continuum in the energy absorption spectrum denotes a 1 — » 3 
p-h continuum and the kink in the middle is indicative of an 
exciton lying inside it. Both the continuum and the kink ap- 
pear in all of the inter-subband transition rate plots in this case. 

In the remainder of this paper, we explore the ideas summa- 
rized above in detail. This paper is organized as follows: the 
microscopic model is introduced in Sec. [IT] and the Hartree- 
Fock equation is discussed in Sec. Ill The numerical results 
obtained by solving the Hartree-Fock equation is discussed in 



Sec. Ill A Lattice modulation spectroscopy experiments are 



reviewed in Sec. IV and the calculation of energy absorption 
spectrum and inter-subband transition rates in the TDHF ap- 



proximation is discussed in Sec. |IV A The numerical results 
are presented and discussed in Sec. IV B The paper is con- 



cluded by a short discussion on the experimental outlook of 
the presented results. 



II. THE MICROSCOPIC MODEL 

In this section, we review the microscopic model for spin- 
less fermions with electric dipole-dipole interactions in a one- 
dimensional optical trap. This model is relevant to a pancake 



of polar molecules prepared in a single hyperfine state, as well 
as to a stack of well-separated pancakes. 

For concreteness, we assume that the gas is confined about 
the x-y using a confining optical potential centered at z = 
0. Also, we assume that the dipoles are aligned along the z- 
axis (perpendicular to the confining plane) using a strong d.c. 
electric field. The Hamiltonian of the system is the sum of 
the optical trap potential, the kinetic energy and the electric 
dipole-dipole interactions. A convenient basis for the second 
quantized notation is one that diagonalizes the one-body part 
of the Hamiltonian. We choose the following basis: 



(r|a,k> 



1 

(t>a{z)—=e l 
VA 



(2) 



where 4> a (z) is the wavefunction of the a'th transverse mode 
of the trap, A is the area of the trap in the x-y plane and x = 
(x, y) is the in-plane coordinates. In this basis, the second- 
quantized Hamiltonian is easily found to be: 



k , a 



h 2 \k\- 

2m 



e « I C k.a C k,< 



+ E V «/5;7A(q)c kl +q, Ct C L-q, 7 C k 2 ,A C k 1 , / 3' ( 3 ) 

a/3; 7 A kik 2 q 

where Ck, Q (c^ ) annihilates (creates) a particle in the in- 
plane momentum state k and a'th subband. e a is the zero- 
point energy of a'th subband (a = 0, 1, . . .). We measure the 
energies with respect to the zero-point energy of the lowest 
subband (a = 0). V a p n \(q) is the Fourier transform of the 
effective inter-subband dipolar interaction, defined as: 



V" q /3; 7 a(x-x') = / / dz dz' 4>* a (z) <f>*(z') <f>\(z') <j)p(z) 

x V dip (z-z',x -*.'), (4) 

where Vdi p (z,x) is the electric dipole-dipole interaction of 
two particles with a center separation of (z, x): 



y di p(z,x) = 



D 2 



(|x| 2 + z 2 ) 



r(l*l : 



2z l 



(5) 



The intra-subband interactions are special cases of Eq. Q, 
i.e. the intra-subband interaction in a'th subband is given by 
V Q a;aa(q)- Since the dipoles are perpendicular to the x-y 
plane, the interaction matrix elements are isotropic in the in- 
plane momentum q. 

One can always choose the transverse wavefunctions 
4> a (z) to be real and of well-defined parity for trap potentials 
which are symmetric about z — (harmonic trap is one 
example). As a result, it is easy to verify that V q /3 ;7 a 
is invariant under the following interchange of indices: 
a o (i, 7 «-» A, or af3 o 7A. Also, V a /3 n \ vanishes if 
a + /3 + 7 + A= 1 (mod 2). This parity-conserving behavior 
is due to the invariance of Eq. |5) under inversion z — > —2, 
and breaks down for unaligned dipoles. 
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For simplicity, we assume that the optical trap is perfectly 
harmonic with a frequency u; trap . In this case, we easily find: 

€ a = a (/iWtrap), 

<j> a {z) = -= 1 H a (z/ ai _)e~ z2 / 2a ^, (6) 

Vtt 1 / 2 a! 2 a a± 

where H a (z) is the a's Hermite polynomial [50 1 and a± is the 
transverse confinement width, related to the trap frequency as 
^trap = h/ma\. A generating function for V a i3 n \(q) for 
harmonic traps is derived in Appendix |A| and explicit expres- 
sions for the first few inter-subband interactions are given. 

We note that in the experiments, nearly-perfect harmonic 
trapping can be achieved by loading the gas into a single well 
of a strong optical lattice. An optical lattice potential of the 
form V\ a t.( z ) = {ti 2 /2ma\k 2 ) sin 2 kz yields a harmonic trap 
with frequency h/ma 2 ^ centered at z = in the limit k — » 0. 



HI. HARTREE-FOCK EQUATIONS FOR A UNIFORM 
QUASI-TWO-DIMENSIONAL FERMIONIC GAS 

We briefly review the Hartree-Fock (HF) theory for a 
uniform quasi-two-dimensional fermionic gas in thermal 
equilibrium and apply the formalism to the system of 
fermionic polar molecules. The Hartree-Fock equation 
for quasi-two-dimensional systems is found to be sub- 
stantially more difficult to solve compared to the strictly 
two-dimensional case due to subband hybridization. 

We begin our treatment with the usual definition of the 1 - 
particle thermal Green's function: 



0h 



dre^ lV[p GCk „(T)q 



(7) 

where /3 — 1/ksT, T is the temperature, icj n = ( 2n + 1 ) 7r//3 
is the fermionic Matsubara frequency, pq = er^ H /Tr[e~^ ff ] 
is the grand-canonical statistical weighting operator and 

„(t) 



(t) = e rH / h e£ e rH l h is the imaginary-time Heisen- 



berg fermion annihilation (creation) operator. In the pres- 
ence of interactions, the non-interacting subband indices no 
longer remain good quantum numbers due to hybridization 
and consequently, (^(k, iu> n ) is expected to have non-zero 
off-diagonal elements. 

The Hartree-Fock approximation for the Green's function is 
given by the following diagrammatic Dyson's equation |51 1: 
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+ k' 



;q-k' (8) 



where the thin and thick fermion lines denote non-interacting 
and interacting Green's functions respectively. The diagram 
yields the following equation: 



g a p(k,iu n ) = S2 1 9(k,tWn) + a^(k > »w n )E* I/ (k) 
x g v pQ&.,iu> n ), 



where the proper self-energy S* (k) is the sum of the direct 
and exchange diagrams: 

KM = \ E / 70^2 [WW - V,A; 7 ,(k - k')] 



(27T) 

x5 Al (k»< 0+ . 



(10) 



Summation over repeated indices is assumed throughout this 
paper. The non-interacting thermal Green's function is given 
by: 



£°^(k, iu n ) 



6 



k.n 



2m 



-e a -fx. (11) 



In the absence of spontaneous symmetry breaking of the inver- 
sion symmetry (which may happen in the strongly interacting 
regime and is beyond the scope of this paper), the interactions 
only mix subbands of the same parity due to the symmetries 
of the inter-subband interaction matrix elements mentioned in 
the previous section. Investigating Eqs. (|9]l and (lOi shows 
that Q a /3 = 0ifa + /3 = 1 (mod 2) and therefore, Q a p will 
have the following matrix structure in the non-interacting sub- 
band indices: 



Q = 



( £oo 




V : 





Gii 


Gl3 





G02 


G22 


G24 





Gl3 



£33 




G04 




£44 



(12) 



In order to solve the Hartree-Fock equation at finite tempera- 
tures, one must carry out the Matsubara summations appear- 
ing in Eq. (10 1. The summation can be easily done once 
the dependence of the Green's function on the Matsubara fre- 
quency is explicitly known. In fact, expressing the interacting 
Green's function its diagonal basis reveals its explicit depen- 
dence on the Matsubara frequency. Treated as a matrix equa- 
tion, Eq. Q formally yields: 



(13) 



G = 



- Q/ 3( k ) = s a/?( k ) + £k, Q <W, i-e. the proper self- 



where E c 

energy including the kinetic energy contribution. It is easy 
to see that E (k) is a symmetric matrix and therefore, there 
exists a real orthonormal basis in which it is diagonal and has 
real eigenvalues. Let U(k) be the unitary transformation that 
diagonalizes E (k): 

E*(k) = I7(k)S(k) t/ T (k), 
S(k) = diag{&(k), | 2 (k), ...}, (14) 

where {£ Q (k)} are the eigenvalues. The same transformation 
clearly diagonalizes the interacting Green's function: 

G a ^,iu n ) = [U(k) T GU(k)] ap = [iw n l - S(k)]^ 

5 a p 



(9) 



(15) 
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As promised, the Matsubara frequency summation appearing 
in Eq. < fT0] > can be evaluated with ease in the new basis: 

= ^EE C7 v(k / )Kr P (k / 



= J2Ux P (V)n F lUk>)}U JP (k>). 



(16) 



where n F (x) = (e l3x + l) 1 is the Fermi occupation func- 
tion, resulting from the summation over fermionic Matsubara 
frequencies BP . Plugging this result into Eq. ([9]), we get an 
explicit self-consistent equation for the proper self-energy: 



S*,(k) = J 



d 2 k' 



DV ;7 a(0) - V MA;7 „(k - k')] Ux p Qs!) 



x U lp {k')n F fo(k') 



(17) 



It is understood that £/(k) and £ p (k) are implicit functions of 
£*(k), defined in Eq. M. 



It is clear from the preceding discussion that the orthogo- 
nal transformation J7(k) also defines the mean-field single- 
particle states in the presence of interactions. We identify 
£ a (k) as the energy dispersion of a'th Hartree-Fock subband 
and define the Hartree-Fock fermion annihilation (creation) 
operators as: 



4!L(T) = E^( k ) c k]>)- 



(18) 



It is straight-forward to show in light of Eq. (15 1 that 
Q pv (k, ioj n ) can be identically defined using the Hartree-Fock 
fermion operators: 



G M v(k,iu n ) = - 



(19) 

where /5q F is the grand-canonical operator defined in terms 
of the Hartree-Fock decoupled Hamiltonian. 

At this point, it is also useful to define the effective inter- 
action between Hartree-Fock quasiparticles. The evaluation 
of response functions, which is our goal in the next section, 
is more natural in this basis. Expressing the interaction part 
of the Hamiltonian in terms of Hartree-Fock fermion opera- 
tors, one can easily read off the renormalized interaction be- 
tween incoming particles in Hartree-Fock subbands f3 and A 
with momenta ki and k2, scattering to subbands a and 7 with 
momenta ki + q and k2 q respectively: 

Va/3; T A(kl,k 2 ,q) = ^ V a '/?' ;7 'A' (q) U a > a {ki + q) 

x L/y 7 (k 2 - q) U x >x(k 2 ) CWkx). (20) 

Note that the renormalized interaction is no longer just a func- 
tion of the momentum transfer, but also depends on the indi- 
vidual momenta of the scattering quasiparticles. 



A. Results: renormalized bands 

In this section, we present the numerical results obtained 
by solving the Hartree-Fock equation. The numerical method 
is described in Appendix [B] in detail. The temperatures are 
reported in the units of T*, which is defined as: 



T* = h 2 n/2mk B . 



Note that T* ua related to the Fermi temperature of a 
two-dimensional non-interacting Fermi gas, as: 



(21) 



T* =rj, 0) / 47r -0.08T 



(0) 



(22) 



Fig. [2] shows the energies of the first five Hartree-Fock 
subbands (panels a, c and e) and their corresponding den- 
sity of states (DOS) (panels b, d and f) for = 1.0, 
T/T* = 0.1 and for three different transverse confinement 
widths ^/na± = 0.2, 0.4 and 0.6. Deviations from the non- 
interacting quadratic energy dispersions can also be observed 
in the DOS plots: the DOS of a quasi-two-dimensional non- 
interacting Fermi gas has a uniform staircase structure (shown 
in Fig. [2J3, and^ as blue dashed lines for reference). In 
the presence of interactions, we find that (1) the DOS plot 
starts at a finite energy, meaning that the zero-point energy of 
the lowest subband is lifted, (2) the DOS plot no longer has 
flat regions (due to deviations from quadratic), (3) the energy 
spacing between the jumps, which correspond to the spacing 
between the zero-point energies of the subbands, become non- 
uniform, and (4) the DOS of a non-interacting gas is always 
larger than the DOS of the interacting gas, which is associ- 
ated to the long-range repulsive nature of dipolar interactions 
in the studied confined geometry. 

At zero temperature, the occupation of higher subbands is 
only due to Pauli exclusion and the Fermi occupation func- 
tion is sharp. Therefore, at any given density, trap strength 
and interaction strength, only a finite number of Hartree-Fock 
subbands are fully or partially occupied. There, a phase dia- 
gram can be obtained for the system at T = as a function 
of trap and dipolar interaction strengths (Fig.|4|. It is noticed 
that at fixed transverse confinement width a± and density n, 
stronger interactions result in occupation of higher subbands. 
This behavior can be understood in light of the reduction of 
DOS due to interaction. The energy and fractional occupation 
of the subbands were found to be continuous across the phase 
boundaries and therefore, the transitions are continuous. The 
same transitions has been reported to be first-order for quasi- 
two-dimensional electron gas [52|. 

In order to see the hybridization of non-interacting sub- 
bands, we have plotted the fractional density of non- 
interacting and Hartree-Fock subbands in Fig. [3^ and [3j? re- 
spectively, as a function of transverse confinement width and 
at fixed r d = 0.05 and temperature T/T* = 0.1. Fig. 
and |3}; show the same quantities at a higher temperature 
T/T* — 1.0. It is observed that larger transverse confine- 
ment widths (i.e. weaker traps) naturally results in occu- 
pation of higher subbands. Hybridization is clearly notice- 
able by comparing Fig. [3^ and[3]3: occupation of the lowest 
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FIG. 2. The energy dispersion of the first five Hartree-Fock subbands (panels a, c, and e) and the density of states (DOS) g(e) (panels b, d and 
f) for r d = 1.0, T/T* = 0.1. (a) and (b): ^na± = 0.2 (c) and (d): ^na± = 0.4 (e) and (f): ^/nax = 0.6. The DOS in the absence of 
interactions is also shown for comparison (blue dashed lines). 
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FIG. 3. Subband mixing and renormalization of the chemical potential. Panels (a) and (b) show the fractional population of non-interacting 
subbands and Hartree-Fock subbands respectively, as a function of transverse confinement width a± at constant dipolar interaction strength 
Td — 0.05 and temperature T = 0.1 T* ~ 0.008 Tp\ Panels (c) and (d) are the same quantities at a higher temperature T — T* ~ 
0.08 t£, ' . Panel (e) shows the chemical potential as a function of transverse confinement width a± at T = 0. 1 T* ~ 0.008 Tp ' (solid blue), 
T = 1.0 T* ~ 0.08 T<. 0) (solid black) and T = 2.0 T* ~ 0.16 T { F 0) (solid red). The dashed lines denote the chemical potential in the absence 
of interactions for comparison. Notice the non-monotonic ordering of the chemical potential vs. temperature in the inset (refer to the text for 
discussions). 



Hartree-Fock subband amounts to occupation of several non- 
interacting subbands of even parity, 0,2,4, . . . , with decreas- 
ing weights. A consequence of subband mixing is anomalous 
population inversion, i.e. the 2nd non-interacting subband is 
populated before the 1st non-interacting subband, due to hy- 
bridization with 0th, 4th, . . . subbands. 

Finally, Fig. [3^ shows the renormalized chemical potential 
as a function of transverse confinement width at three dif- 
ferent temperatures (solid lines). We have also plotted the 
chemical potential in the absence of interactions for reference 
(dashed lines). It is noticed that the renormalized chemical po- 
tential is always greater than its non-interacting value, which 



is again due to reduction of the DOS in the presence of in- 
teractions. An interesting observation is the non-monotonic 
behavior of the chemical potential as a function of tempera- 
ture, which is clearly noticeable in the inset plot of Fig. [3^. 
As a consequence, the isothermal compressibility of the in- 
teracting gas, kt = n~ 2 (dn/dn)T also turns out to have a 
non-monotonic behavior as a function of temperature. This 
behavior has been reported earlier [53] and is confirmed by 
our calculations. Intuitively, a small rise in temperature will 
result in thermal excitation of states just below the Fermi 
level. Having larger momenta, the thermally excited quasi- 
particles experience stronger interactions and decrease DOS 
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FIG. 4. The quasi-two-dimensional liquid phases of fermionic polar 
molecules at T — 0. The abbreviations No, Ni, etc stand for nor- 
mal liquid phases populating up to the zeroth, first, ... Hartree-Fock 
subbands respectively. 



at the Fermi level. Consequently, the chemical potential has 
to be increased in order to compensate for the reduced DOS. 
At higher temperatures, the effect is suppressed as the states 
above the Fermi level are significantly populated and a smaller 
chemical potential is required to keep the number of particles 
constant. 



IV. PROBING THE INTER-SUBBAND EXCITATIONS IN 
LATTICE MODULATION SPECTROSCOPY 
EXPERIMENTS 

The single particle and collective excitations of a strictly 
two-dimensional gas of polar fermionic molecules has been 
recently studied by several authors 11221 l33l [361 . The two- 
dimensional gas corresponds to the single-subband limit of a 
quasi-two-dimensional gas and can be achieved by increasing 
the trap frequency. At zero tilt angle, i.e. when the dipoles 
are aligned perpendicular to the two-dimensional plane of the 
trap, the dipole-dipole interactions are effectively repulsive. 
The dipolar interactions are also short-range and regular in 
the long wavelength limit. Therefore, such systems are ex- 
pected have the same qualitative properties as He-3, which is 
described well by Fermi liquid theory of neutral systems ll54l . 
The elementary excitations of such systems consists of single 
particle excitations and the zero sound. 

As the higher subbands are populated, we encounter a new 
class of elementary excitations. Analogous to the strictly two- 
dimensional case, the inter-subband excitations also come in 
two flavors: inter-subband quasiparticle-like excitations and 
p-h bound states (excitons). 

The inter-subband excitations can be experimentally de- 
tected using the technique of lattice modulation spec- 
troscopy 1 44 -48 1 which was originally introduced and real- 
ized in cold atomic systems in order to study the Mott insu- 
lator and superfluid phases of the simulated Hubbard model. 
This method relies on the high controllability of the optical 
lattices which allows rapid modulation of the optical poten- 
tial. Introducing a weak oscillatory modulation to the ampli- 
tude of the optical potential, one stimulates the transition of 
particles from lower subbands to higher subbands. The en- 
ergy of the inter-subband excitations must match the modu- 



lation frequency for the transition to occur. Therefore, one 
directly measures the excitation energy in such experiments 
and no tedious calibration of parameters is needed. 

If one is interested to measure the energy absorbed in the 
process of lattice modulation, one allows the system to re- 
thermalize after the lattice modulation pulse. All optical po- 
tentials (including the polarizing d.c. electric field in the case 
of polar molecules) are suddenly switched off then and the 
density profile of the gas is measured after an interval of bal- 
listic expansion. The change in the width of the central peak 
in the momentum profile can be taken as a measure of the 
absorbed energy. By carrying out this procedure for a range 
of frequencies, one obtains the inter-subband excitation spec- 
trum of the system. We refer to this procedure as energy- 
resolved measurement for brevity. 

Another quantity which is often measured in lattice modu- 
lation spectroscopy experiments is the change in the number 
of particles in each subband. We refer to such experiments as 
band-resolved measurements for brevity. In such experiments, 
one avoids the re-thermalization of the gas after the modula- 
tion pulse. Instead, the gas is allowed to expand ballistically 
quickly after the modulation pulse and the transverse momen- 
tum profile of the particles is measured. The particles are re- 
solved into subbands by fitting the measured density profile 
to a weighted sum of the density profiles calculated from the 
transverse wavefunctions. Band-resolved measurements can 
be used to study the mechanism of energy absorption. 



A. The Perturbative Formulation of Lattice Modulation 
Spectroscopy Experiments 

In this section, we describe a theoretical framework for pre- 
dicting the results of modulation spectroscopy experiments of 
polar molecules. As mentioned in Sec. [TTJ we are interested 
in the case of well separated layers, where one can neglect 
the inter-layer couplings and just focus on transitions within 
a single layer. For concreteness, we focus on the pancake at 
z = 0. Upon introducing the amplitude modulation, the opti- 
cal potential takes the following time-dependent form: 



Mat.M) 



Wo 

1 



2 2 
2 mW trap^ 



Vtfit)] sin 2 (fcz) 



(23) 



where w trap = ky / 2V /m, a±_ = y/h/(mu} tra . p ), £ = 
(fiwtr a p/2Vb)Vi and f(t) is the shape of the amplitude mod- 
ulation pulse. We have expanded the optical potential to 
quadratic order about its minimum at z = 0. The second 
quantized form of the lattice modulation potential, i.e. the 



second term in Eq. ( 23 1, can be written as 



*i(t) = £/(*) _ 

k,a/3 



^ Te»/3 c ka c k,8> 



(24) 
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where T 



a/3 



j dz 4>* a {z) z 2 4>p{z). For a harmonic trap, 



T Q) 3 is given by: 

T Q(3 = I( A /( a + l)(a + 2) 5 a ^ 2 

+ y/a(a - 1) 5 a . p+2 + (2a + 1) 5 a p] . (25) 

It is easier for subsequent derivations to express the perturba- 
tion Hamiltonian in terms of Hartree-Fock fermion operators: 



where: 



T Qj a(k) n Q /3(k) 

k,a/3 



J-a/3 a'/3' 
a'/3' 



(26) 



(27) 



We note that both T Qj g and T a p are symmetric is the 



subband indices. It is clear from Eqs. (24i and ( |25] l that 
this perturbation stimulates an inter-subband p-h excitation 
between subbands whose index differ by ±2. Note that we 
have tacitly assumed that the amplitude of the modulation 
is low (i.e. £ <C 7ik>trap in Eq. |23] > and neglected the small 
anharmonic terms such as ~ cos(ut) z 4 , which stimulate 
transitions between subbands whose index differ by ±4 as 
well. 

For low-amplitude modulations and short modulation 
pulse durations, the disturbance in the equilibrium state of 
the system is expected to be negligible and therefore, we 
can treat the modulation potential in perturbation theory. In 
the first order approximation, the response of any dynamical 
variable of system to an oscillatory external field is also 
purely oscillatory. Therefore, the linear response vanishes 
on average and one needs to consider processes which 
are at least second order in the external fields in order to 
describe transport phenomena such as energy absorption and 
inter-band transitions. 

We focus on energy-resolved measurements first. A 
straightforward calculation yields the following expression for 
E(t) = dtTr[p(t) H] in the second-order perturbation theory: 

E(t) = 2<Ke j At' (V moA (t)HV mod {t r ) 

J—oot 

-HV mod (t)V mod (t')) 



iff(t) dt'f(t')d t T(t-t'), 



where: 



and: 



(28) 



E E W(ki,k2,f), (29) 

ki ,CK/3 k2 ,"f\ 



f a/ 3(ki) f 7A (k 2 ) n^. 7A (k l5 k 2 ; t), 
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FIG. 5. Plot of the a.c. pulse for rj/fl 
measurement is done at t = 0. 



0.05. The band-resolved 



in which IT^ /3 . 7A (ki, k 2 ; t) is the retarded inter-subband po- 
larization insertion: 

ri^j. 7A (ki , k 2 ; t) = -i9(t) Tr { Po [n 0a (ki , t) , n lX (k 2 , 0)] } . 

(31) 

In the experiments, the lattice modulation pulse is switched 
on and off smoothly. For concreteness, we assume an expo- 
nentially switched-on pulse: 



/(i) = e-"l*l cos(m) 



(32) 



where 77 is the switching rate and £1 is the frequency of the a.c. 
modulation. We chose the exponentially switched-on pulse 
(shown in Fig. [5]l since it yields simple and transparent ana- 
lytical expressions. It is shown in Appendix |C] that the energy 
absorption rate for such a pulse is given by: 



E\ n n = ~-flJm [T(n + ir,)] +0(he r^-i p ). 



(33) 



Note that in light of Eqs. (29 > and (30 1, T(w)> the Fourier 



transform of Tit), can be identified as the self-energy cor- 
rection of the optical lattice photons with energy Hu> coupled 
to the fermionic molecules. Thus, Eq. (33i can be simply 



interpreted as the energy of the absorbed photons multiplied 
by their decay rate, i.e. the imaginary part of their self-energy. 



It is evident from Eq. ( |29[ > that the fundamental quan- 
tity to be evaluated is A (ki,k 2 ;t), i.e. the dynami- 
cal inter-subband polarization insertion with zero net momen- 
tum transfer. We evaluate this quantity in the conserving 
time-dependent Hartree-Fock (TDHF) approximation (49), 
also known as Generalized Random Phase Approximation 
(GRPA) ||54l or Random Phase Approximation with Ex- 
change (RPAE) |55|. Once the inter-subband polarizations 
are found, one can easily evaluate the energy absorption rate 
using Eqs. ( f2"9"l >, p0| ) and ( [33| . We carry out the diagram- 
matic calculations in the imaginary-time formalism. The re- 



tarded polarization appearing in Eq. ( 30 1 can be found by the 
standard procedure of analytical continuation of the Matsub- 
ara frequency to the upper complex frequency half-plane: 



n a /3;7A(ki,k 2 ;f2) = n Q/3;7A (ki,k 2 ; iv n -> SI). 



(34) 



(30) 



Evaluating polarization diagrams in the TDHF approxima- 
tion amounts to summing the direct and exchange scatterings 
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between the p-h pairs to all orders in each p-h loop. In other 
words, the full polarization diagram is the sum of all ring di- 
agrams with ladder-like vertex corrections |j49l . The succes- 
sive iterations of following Bethe-Salpeter equation generates 
all such contributions: 



ki,/3 



k 2 ,\ 



k\.J5 k' .v 



k 2 ,X 



6k 



ki fc 2 



(35) 



ki,a 



k 2 ,l 



h\,ot fei,a k' ,fi 



k 2 ,l 



where the double zigzag line represents the combined direct 
and exchange interaction between Hartree-Fock quasiparti- 
cles: 



k , a 



k',v k,0 



k' ,/j k,a 



k'v 



k' ,/i 



k,/3 k',v 
JVWWi' 

..>'' v.. 

k,a k' ,[i 



= V^. a/i (k, k', k' - k) - V a ^(k, k', 0). 

(36) 

The negative sign of the direct interaction is due to the extra 
Fermion loop it introduces. The diagrammatic Eq. ( 35 1 yields: 

lj2 r d2k ' 



i j ( 27r )2^STi( k i> k ') 

x fl^ n \(k', k 2 ; iv n ) 



(37) 



where k\ — (ki, iui ni ), k 2 = (k2,iw ri2 ) and k! = (k',iuj' n ). 
We note that if we had formulated the problem in the non- 
interacting basis, the Green's functions would have been non- 
diagonal in the subband indices and the resulting integral 
equation would have had additional intermediate subband in- 
dex summations. 

Summing both sides of Eq. ( [37] i over the Matsubara fre- 
quencies ico ni and iuj n2 an d analytically continuing iv n — > Q, 
we get: 



n a/3;7A (k 1 ,k 2 ;0)=fl^(k 1 ;0) 



d 2 k ; .-(d+e 
(2tt) 2 a ^ 



CT^(ki,k')n^ ;7A (k',k 2 ;Q) 



(38) 



where, nj^k; f2), the bare inter-subband polarization with 
zero net momentum transfer, is defined as: 



np(k) - n a (k) 

n + ^(k)-| a (k)* 



(39) 



The poles of n^2(k; fl) correspond p-h excitation 



energies. 



Eq. ( 38 i is an integral equation for the inter-subband polar- 



ization. However, it is more suitable for our purpose to formu- 
late an integral equation for the following auxiliary function, 



in which some of the summations appearing in Eq. ( 29 1 is al- 
ready carried out: 

r „ Ek I „Ai>ft,A(ki,k 2 ;S!)T, A (k 2 ) 

"° 0<k,;n, = 5>^> ' <40> 



Using Eq. 1(38} , one easily finds: 

ftojsCkijfi) =T a /s(ki) 
d 2 k' 



, 2;7)2 «$$(ki, k') fiW (k' ; 0) V(k'; O). 



(41) 



Iterating the above integral equation, 




one identifies h a p as the effective inter-subband transition 
matrix. The presence of p-h scatterings to all orders suggests 
that ft Qj a(k;f2) will have poles at excitonic energies. Note 
that in contrast to Eq. ( 38 1, the bare polarizations (which have 
integrable singularities) appear under integrals in the above 
equation. Therefore, provided that its Fredholm determinant 
is non-vanishing, Eq. ( [41} yields smooth solutions. However, 
the Fredholm determinant may vanish for an isolated set 
of frequencies, corresponding to the excitons. We use this 
criteria for locating the excitons. 

Once ft, Q fl(k;Q) is found, T(f2) can be readily evaluated 
using Eq. ( ftO} : 

T(fi) =2Z)T a/ ,(k 1 )nW(ki;n)^(k 1 ;n). (42) 

a/3 ki 

The numerical procedure used to solve the integral equation 
for h a p is described in detail in Appendix^ 

In the TDHF approximation, the excitons that lie outside of 
p-h continuums are undamped and have infinitely long life- 
times. Therefore, the excitonic poles appear infinitesimally 



below the real frequency axis. Thus, according to Eq. (33i 



such undamped excitons yield sharp Dirac delta-like peaks in 
the energy absorption spectrum in the adiabatic pulse switch- 
ing limit (7/ — > 0). The spectral weight of an exciton can 
be found using the analyticity of T(w) in the upper complex 
frequency half-plane. Assuming that that T(w) is singular at 
uj = ljq — i0 + , one can find the spectral weight associated to 
it using the Kramers-Kronig relations: 



W 



4 



uj lim I (w - ui ) 9\e\T(uj)\ }. 

oj—mj I ) 



(43) 
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Finally, we address the subtler problem of evaluating the 
inter-subband transition rates. The subtlety originates from 
the fact that the number of particles in each subbands is not 
conserved separately by the microscopic Hamiltonian. The 
presence of inter-subband interaction matrix elements such as 
V 7Q ;7/3 implies that [N y ,H] ^ 0, where N y is the total num- 
ber operator of particles in 7'th subband. We note that this 
result holds for Hartree-Fock number operators as well. As 
a consequence, the inter-subband excitations generated by the 
lattice modulation term are subject to non-equilibrium dynam- 
ics as the system evolves. Ultimately, one needs to use the 
Schwinger-Keldysh formalism and non-equilibrium Green's 
functions in order to calculate the transition rates in such sit- 
uations. We note that this additional difficulty is not present 
in theories in which the number of each species is conserved 
separately. In those cases, the transition rates can be evaluated 
using the equilibrium formalism in the second-order perturba- 
tion theory |56|. 

We quote the final result of the non-equilibrium analysis 
and leave the technical details for future works. In brief, our 
strategy is to solve the non-equilibrium Dyson's equation up 
to second order in the external field. The inter-subband transi- 
tion rates can be extracted from the second-order corrections 
to the non-equilibrium Green's functions. For the exponen- 
tially switched-on pulse, we find that the total number of par- 
ticles transferred from the Hartree-Fock subband a to /? for 
times t < can be written as: 



,—2771*1 . 



(<0 



\hp a (n+iTi)\'3m U { °hn + ir]) 



(44) 

We have neglected the oscillatory terms in the above equa- 
tion. Such contributions are the second harmonics of the ex- 
ternal a.c. field which naturally arise in the second-order per- 
turbation theory. In the limit rj/Q, -C 1, such terms are much 
smaller than the retained term. Moreover, they yield no steady 
transition current on average due to their oscillatory nature. 

The effective duration of the exponentially switched-on 
pulse from t — — 00 to t = 0, when the measurement takes 
place, is At ~ -q^ 1 . We define the effective inter-subband 
transition rate as: 



(45) 



For weakly interacting systems, subband hybridization is 
a negligible effect (see Fig. [3j the fractional hybridization 
is ~ 10 -6 for rd — 0.05) and the transition rates be- 
tween Hartree-Fock subbands are virtually the same as the 
transition rates between the non-interacting subbands, i.e 

We conclude this section by a short discussion on the condi- 
tions for validity of the time-dependent Hartree-Fock theory. 
We begin with the subtler case of evaluating the rate of change 
of non-conserved quantities. Note that conserved and non- 
conserved quantities are defined with respect to the Hamilto- 
nian in the absence of external fields. Generally, we expect the 



prediction of the Hartree-Fock theory to become less reliable 
for longer pulse durations. In particular, the study of the adia- 
batic limit must be avoided in this approximation. The reason 
is that the collision integrals are neglected in the Hartree-Fock 
approximation and consequently, the resulting picture lacks 
relaxation and rethermalization mechanisms. Hence, the the- 
ory yields unphysically long-lived excitations (with the only 
exception of bound states lying inside p-h continuums) and 
leads to erroneous predictions for transition rates. Therefore, 
one must make sure that the lifetime of excitations are larger 
than the pulse duration. For the exponentially switched-on 



pulse of Eq. ( 32 1, this condition leads to the following crite- 
ria: 



0T-}<2r? (for all a), 



(46) 



where (r^ cl ) is the average relaxation rate of an excitation in 
a'th subband. An estimate of (r^ 1 ) for quasiparticle-like ex- 
citations is provided in Appendix [D] We will use Eq. (46 1 and 
the results of Appendix [D] in order to choose an appropriate 
value for 77 in the next section. We note that in experiments 
with reactive polar molecules, another constraint on the pulse 
duration results from the reduced lifetime of molecules in the 
higher subbands. In practice, 77 must be chosen to be larger 
than both the relaxation rate and the molecule loss rate. 

Finally, we remark that the reliability the Hartree-Fock pre- 
dictions for the evolution of conserved quantities, such as the 
expectation value of energy, is expected to be unaffected by 
long pulse durations. The reason is that the rate of energy ab- 
sorption processes, i.e. creation of excitations, is essentially 
determined by the equilibrium states of the system and their 
occupation, which is assumed to be only slightly disturbed 
by the external perturbation during the experiment. In other 
words, relaxation processes only affect the evolution of non- 
conserved quantities. 



B. Results 



In this section, we present the results obtained by evaluating 



Eq. ( 33 1 and ( 44 1 numerically. We keep the first five subbands 



in the numerical calculations. The results are given in terms 
of the following dimensionless and intensive quantities: 



i= {2hN- 1 C 2 )E, 



w 



(2hN~ 1 C 



trap 



)W, 



fi = fJ/Wtrap- 



(47) 



Based on the remarks mentioned at the end of the previous 
section, one expects the energy absorption spectrum obtained 
in experiments done at a finite switching rate to be essentially 
a broadened version of those obtained in the adiabatic limit. In 
fact, the appearance of T(f2 + irf) in Eq. (33 1 and the analyt- 
icity of T(w) in the upper complex frequency half-plane is a 
rigorous justification of this claim. In light of this observation, 
we have evaluated the energy absorption rates in the adiabatic 
limit for the clarity of presentation. In this limit, the excitonic 
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peaks will be Dirac delta-like and one can easily differentiate 
between the excitons and the p-h continuums. 



1. Energy Absorption Rates 




FIG. 6. Energy absorption rate plots for various dipolar interaction 
strengths at fixed ^na± = 0.2andT = 0.1T* ~ 0.008 T ( F 0) . In all 
of the shown cases, the spectrum consists of a single p-h excitation 
continuum (yellow shaded continuum) and a single anti-bound exci- 
ton (red spikes). The exciton captures 99.74% of the spectral weight 
for the range of rd shown in the figure. 

We start with the simpler case of a strong trap at low tem- 
peratures, where only the zeroth subband is populated. Fig. [6] 
shows the energy absorption spectrum for a range of weak 
dipolar interactions r d = 0.001, 0.005, 0.010, . . . , 0.050 
at a constant transverse confinement width y/na± = 0.2 and 
temperature T = 0.1 T* ~ 0.008 t| 0) . It is noticed that the 
spectrum consists of a single p-h continuum, corresponding 
to — > 2 transitions, and a single exciton. The energy of the 
exciton lies above the p-h continuum and therefore, it is an 
anti-bound p-h pair. The anti-binding nature of the excitonic 
mode is a consequence of the anisotropic dipole-dipole 
interactions. The exciton captures 99.74% ± 0.01% of 
the spectral weight, leaving only 0.26% ± 0.01% for the 
continuum. These ratios were found to be independent of rd 
in the studied range, though, may depend on a± and T. 

Fig. [7] shows the energy absorption spectrum for various 
temperatures in the range 0.04 < T/Tp^ < 1.60 at a con- 
stant transverse confinement width y/na± = 0.2 and dipolar 
interaction strength = 0.05. It is found that at low temper- 
atures, the spectrum consists of a single exciton and a single 
p-h continuum corresponding to — > 2 transitions (yellow 
regions). As the temperature is increased, the spectral weight 
of the excitonic mode slightly drops from 99.74% ± 0.01% to 
99.67% ± 0.01%. At this point, the p-h continuum associated 
to 1 — » 3 transitions (shown as blue) becomes visible due 
to thermal population of the first excited subband. The 
exciton also merges into the 1 — s- 3 continuum. At higher 
temperatures, more p-h continuums appear continuously due 




FIG. 7. Energy absorption rate plots for various temperatures at fixed 
\/na± = 0.2 and r<j = 0.05 (achievable in the current experiments 
of the group at IILA). The excitons appear as red spikes and their 
fractional spectral weights are shown above them. At low tempera- 
tures, the spectrum consists of a single p-h continuum corresponding 
to — > 2 transitions (yellow shaded continuums) and a single exciton 
(red spikes). At higher temperatures, more p-h excitation continuums 
appear continuously due to thermal occupation of higher subbands. 
The blue and violet shaded continuums correspond to 1 — > 3 and 
2 — > 4 transitions respectively. The excitonic modes capture more 
than 99% of the spectral weight, whether they are outside or inside 
p-h continuums. 



to thermal population of higher subbands. It is also noticed 
that the 1—^3 exciton, which appears initially on the edge of 
1 — > 3 continuum (the peak at the end of 1 — > 3 continuum, 
visible for T/Tp ~ 0.4 to ~ 1.0), moves out of continuum 

and merges into 2 — !• 4 continuum (visible for T/Tp^ > 1.0). 
We expect that all of the excitons and continuums converge 
to a single peak at u> = 2a; trap at temperatures well beyond 
quantum degeneracy, where all interaction effects are masked 
due to thermal fluctuations. The shift of —> 2 continuum 
towards u> — 2a; tr a P and its shrinking for T > T^ ' agrees 
with this speculation. 

Fig. [8] shows the energy absorption spectrum for various 
transverse confinement widths at a fixed dipolar interaction 
strength r d = 0.05 for two different temperatures T — 
0.1 T* ~ 0.008 r| 0) (a) and T = 10 T* ~ 0.8 (b). The 
low temperature and high temperature results have strikingly 
similar features. At strong confinements (small a±), the spec- 
trum consists of a single exciton and a single p-h continuum. 
Like before, the exciton captures more than 99% of the spec- 
tral weight. Upon relaxing the trap (increasing a±), the energy 
gap between the subbands is reduced and the higher subbands 
will be populated, resulting in appearance of more p-h contin- 
uums and excitonic modes associated to each continuum. In 
the low temperature case (Fig. [8^), the excitons lie outside of 
continuums in most of the plots while at higher temperatures, 
the continuums are broadened and the excitons lie inside the 
continuums in most cases. Had we included higher subbands 
in the calculations, the bare exciton appearing in Fig. [HJ) for 
^/na^_ > 0.3 would have appeared inside the 3—7-5 p-h con- 
tinuum. 
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FIG. 8. Energy absorption rate plots for various transverse confinement widths at fixed ra = 0.05. (a) T = 0.1 T* ~ 0.008 T F , (b) 
T = 10 T* ~ 0.8 Tp°\ The excitons appear as red spikes. For strong transverse confinements (small a±), the spectrum consists of a single 
p-h continuum corresponding to — > 2 transitions (yellow shaded continuums) and a single exciton (red spikes). For weaker transverse 
confinements (larger a±), more p-h excitation continuums appear continuously due to the reduced energy gap between the subbands, corre- 
sponding to transitions 1 — > 3 (blue), 2—^4 (violet), etc. The length of the spikes indicate their dimensionless spectral weight, w. The 
excitonic modes capture more than 99% of the spectral weight, whether they are outside or inside p-h continuums. 



We find that as soon as a new exciton appears (upon in- 
creasing ax), it continuously captures the spectral weight of 
the exciton below it as the trap is relaxed further. Thus, one 
may question the nature of the excitons, i.e. whether they 
are associated to certain inter-subband transitions or have a 
mixed nature and inherit the properties of the previous exci- 
tons which were depleted of spectral weight. We address this 
question in the next section where we discuss inter-subband 
transition rates. 

Another interesting finding is the appearance of Fano line- 
shape [57 | in parts of the spectrum. This phenomenon is most 
easily noticeable in the 1—^3 continuums for y/na± > 0.34 
in Fig.|8]3. The missing spectral weight inside the continuum 
is associated to the interference of the p-h excitations and the 
excitonic mode just above it. 

It is worthy of mention that despite the fact that all of 
the results presented so far belong to the weakly interacting 
regime (r^ < 0.05), one finds that the presence of interactions 
dramatically modifies the inter-subband excitation spectrum. 
In the absence of interactions, all of the structures shown in 
Figs.|6]|7|and|8|disappear, leaving a single peak at u = 2wtrap 
(see FigTj^J. 



2. Inter-subband Transition Rates 

Motivated by band-resolved spectroscopy experiments and 
also in order to elucidate the nature of excitonic peaks in the 
energy absorption spectrum, we evaluate the inter-subband 
transition rates in this section. According to the remarks at the 



end of Sec. IV A and estimates obtained in Appendix |D| for 
the quasiparticle relaxation rate, we choose r\ = 0.01 hn/2m. 



According to Eq. ( 46 1, this choice guarantees that the effec- 
tive duration of the modulation interval is smaller than the 
lifetime of the majority of the inter-subband excitations. 



Fig.|9]shows the energy absorption rate along with different 
inter-subband transition rates for three transverse confinement 
widths at T = OAT* ~ 0.008 T { ° } . The topmost plots are 
slices of Fig. [8^, which are broadened due to finite switching- 
on rate. The energy absorption rates in the adiabatic switching 
limit are also shown for comparison. Note that we have cho- 
sen different scales along the vertical axis for the clarity of 
presentation. 

For y/na± = 0.30 (plots al and a2), the spectrum con- 
sists of a single p-h continuum and an exciton. The only 
non-vanishing inter-subband transition channel is — ^2. For 
y/na± = 0.45 (plots bl, b2 and a3), the spectrum consists 
of two p-h continuums and two excitons. We find that the 
excitons have a mixed nature in this case, i.e. although the 
first and the second excitons are mainly composed of — » 2 
and 1 — > 3 p-h pairs, each has a considerable contribution 
from the other p-h pairs. For y/na± = 0.54 (plots cl, 
c2 and c3), one finds that the spectral weight of the first 
exciton is reduced and second exciton dominates. The en- 
ergy absorption is mainly due to 1 — > 3 transitions in this case. 



Fig. 10 shows the same quantities as Fig. [9] at a higher 



temperature T = 10 T* ~ 0.8T^ 0) . For y^ia ± = 0.05 
(plots al and a2), the spectrum consists of a single p-h 
continuum and an exciton, both of which are associated to 
0^2 transitions. For ^/na± — 0.16 (plots bl, b2 and 
a3), the spectrum consists of two p-h continuums and one 
visible exciton which lies inside the second continuum. The 
exciton is a mixture of both — > 2 and 1 — > 3 p-h pairs. For 
^/na± = 0.28 (plots cl-c4), one finds more continuums and 
excitons. Again, the excitons are mixtures of different p-h 
pairs. 

Another interesting finding is the significantly reduced 
spectral weight between the two excitons in the r ^ 2 plots 
(See Figs. [9]-b2, [9]-c2 and 10 c2). This phenomenon is 
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FIG. 9. Energy absorption rates and inter-subband transition rates for various transverse confinement widths at a fixed dipolar interaction 
strength = 0.05 and temperature T = 0.1 T* ~ 0.008 . The grey shaded plots correspond to quantities evaluated at a finite pulse 
switching rate of -q = 0.01 hn/2m. The blue plots are evaluated in the adiabatic switching limit and are shown for comparison. The red spikes 
correspond to excitons lying outside of p-h continuums in the adiabatic limit. The length of the spikes indicate their dimensionless spectral 
weight, w. (al-2) ^/n~a ± = 0.30, (bl-3) ^na ± = 0.45, (cl-3) ■ s /na 1 _ = 0.58. 




FIG. 10. Energy absorption rates and inter-subband transition rates for various transverse confinement widths at a fixed dipolar interaction 
strength rd = 0.05 and temperature T = 10 T* ~ 0.8 Tp. The grey shaded plots correspond to quantities evaluated at a finite pulse 
switching rate of r\ = 0.01 Tin /2m. The blue plots are evaluated in the adiabatic switching limit and are shown for comparison. The red spikes 
correspond to excitons lying outside of p-h continuums in the adiabatic limit. The length of the spikes indicate their dimensionless spectral 
weight, w. (al-2) ^na± = 0.05, (bl-3) y/na±_ = 0.13, (cl-3) ^/na± = 0.28. 



an example of Fano interference between the excitons and V. EXPERIMENTAL OUTLOOK 

has also been previously reported in semiconductor quantum 

wells f58). Like numerous other theoretical predictions about many- 

body phenomena using polar molecules, the observation of 
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the results presented in this study may also be experimentally 
challenging. In this section, we discuss some of these 
challenges and point out the experimental signatures that are 
expected to more robust. 

The group at JILA have successfully realized a near- 
degenerate gas of KRb polar molecules with a density of 
n ~ 4 x 10 11 cm~ 3 , a temperature of T ~ 1.4 Tp and a 
maximum dipole moment of D ~ 0.22 Debye 151 4151 . Load- 
ing the gas into a one-dimensional optical lattice with a wave- 
length of ~ 1 /im yields stacks of pancakes with y^na± ~ 0.2 
and rd — 0.04, closely matching the parameters used to obtain 
the higher temperature plots in Fig. [7] 

In the experiments with KRb and generally all species 
with energetically favorable two-body chemical reactions, the 
molecules have a finite lifetime due to head-to-tail collisions. 
Populating higher subbands generally results in higher reac- 
tive collision rates, making experiments more challenging. 
From this perspective, it is favorable to prepare the system 
in the single subband regime. The results shown in Fig. [6] 
belong to such regime and in particular, observation of the the 
excitonic peak is expected to be feasible. However, presence 
of heating noise and broadening of excitonic peaks due to 
inelastic scatterings and finite pulse durations may mask the 
p-h continuum, which has less than 1% of the spectral weight. 
If only a single peak is resolved, it can be difficult to judge if 
the peak is associated to the exciton. Therefore, one may need 
to look for signatures of excitons in the multi-subband regime. 

Quemener et al. have recently studied chemical reaction 
rates of reactive polar molecules in pancake geometries and 
have established that the suppression of chemical reactions 
remains effective even if the first few excited subbands are 
populated 11431 . Realizing ultracold gases of non-reactive po- 
lar molecules is another possibility for making quasi-two- 
dimensional systems with long lifetimes. 

In the multi-subband regime, a more robust experimental 
signature for excitonic effects can be found in the plots shown 
in Fig. [8] Appearance of two excitonic peaks (one mainly as- 
sociated to — > 2 and the other to 1 — > 3) and disappearance 
of the first exciton upon relaxing the trap is expected to be 
detectable in both energy-resolved and band-resolved lattice 
modulation spectroscopy experiments. Most strikingly, the 
Fano interference effects lowers the spectral weight between 



the two excitons significantly (see Fig. 10 c2) and allows the 



double peak structure to be resolved even if the broadening 
of peaks is beyond their frequency separation. The plots in 
Fig.JSJ) and Fig. [10] suggest that observation of excitonic ef- 
fects is not limited to strongly degenerate systems and is ex- 
pected to be feasible for temperatures of the order of the Fermi 
temperature as well. 

We also note that extension of the band-resolved mea- 
surement technique to momentum-resolved measure- 
ments ||59l l60ll can be used to directly measure the exciton 
wavefunctions. 

Although the band-resolved measurements yield valuable 
information about the nature of excitations, the energy- 



resolved measurements are expected to be less prone to broad- 
ening since they are not affected by quasiparticle relaxation 
processes. Therefore, energy-resolved measurements are ex- 
pected to yield shaper peaks in the spectrum. 

We also note that we chose the exponentially switched-on 
pulse shape in this study mainly for theoretical convenience. 
Generally, the broadening of the peaks can be minimized by 
ramping up the amplitude of modulations as fast as possible, 
followed by a hold interval at the maximum amplitude, rather 
than a gradual amplitude build-up. 



VI. CONCLUSION 

In this paper, we theoretically analyzed a quasi-two- 
dimensional system of fermionic polar molecules in a har- 
monic transverse confining potential. The electric dipole mo- 
ments of the molecules were assumed to be aligned perpen- 
dicular to the confining plane by applying a strong d.c. elec- 
tric field. We studied the renormalization of the energy bands 
in the Hartree-Fock approximation for various trap strengths, 
dipolar interaction strengths and temperatures (Fig. |2j. The 
renormalized subbands were found to be parity-conserving 
mixtures of the non-interacting subbands (Fig. [3]). A phase di- 
agram was obtained for the normal liquid phases of the system 
at zero temperature (Fig. [4]) as a function of trap and interac- 
tion strengths. 

We also studied the inter-subband excitation spectrum of 
the system in the conserving time-dependent Hartree-Fock 
(TDHF) approximation for various system parameters and 
presented theoretical predictions for energy absorption rates 
in lattice modulation spectroscopy experiments. We found 
that the excitation spectrum consists of both inter-subband p- 
h continuums and anti -bound exciton (Figs. [6] [7] and [8]). The 
excitons capture more than 99% of the spectral weight. Our 
results indicate that both types of excitations are present for 
interaction strengths and temperatures accessible in current 
experiments with polar molecules. 

By evaluating the inter-subband transition rates in lattice 
modulation spectroscopy experiments, we studied the nature 
of excitons and found that they are generally mixtures of 
p-h excitations arising from several subbands (Figs. [9] and 
10 1. We also found a criteria for the validity of the predic- 
tions of time-dependent Hartree-Fock approximation for inter- 



subband transitions (Eq. 46 1 



Finally, we briefly discussed the experimental outlook of 
this study and pointed out the most robust experimental sig- 
natures of excitonic effects based on the the presented results. 
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Appendix A: The generating function for the effective 
inter-subband dipolar interactions in harmonic traps 

In this appendix, we derive a generating function for the ef- 
fective inter-subband dipolar interactions for harmonic traps 
and provide explicit formulas for the first few. We start the 



derivation by transforming the in-plane coordinates of the 
electric dipole-dipole interaction (Eq. [5} to the momentum 
space: 

Vdip(^,q)= /d 2 xe- ic i- x V dip .(z,x) 

( 5(z)-27T J D 2 |q|e~ |q||z| . (Al) 



8ttD< 



The effective inter-subband interaction is obtained by integrat- 
ing out the z coordinate: 



V Q( g ;7 A(q) = / dzdz' 4> a (z) <j>f)[z) (j)j(z') (j)x(z')V dip (z - z' , q) 
8ttD 2 <-°° 



GO pOO 



d?0 a (z)^(2r)0 7 (z')^(z')-47r£> 2 |q| / d£ / 



Mv + O Mv + O 



-2|q|€ 



(A2) 



where the transverse wavefunctions, (j) a (z), etc., are the well- 
known harmonic oscillator wavefunctions (Eq. |6j. We have 
changed variables to r\ — (z + z')/2 and £, = (z — z')/2 in 
order to take care of the absolute value appearing in Eq. < |A 1 [ > 
conveniently. 

A generating function for V Q/ 3 ;7 A(q) can be found by ex- 
pressing the Hermite functions appearing in Eq. (|6|l in terms 



of their generating function: 



E 



H n (x)t n 



(A3) 



Plugging Eq. |6]) into Eq. ( |A2[ ) and using Eq. ( |A3| >, we get the 
following generating function for V a p- n \(q): 
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and V ct(3;7A (q) is given by: 



Using Eq. ( A5 I, we provide explicit expressions for the inter- 



actions between the first two subbands and their long wave- 



V a/9;7 A(q) 
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length limits for reference: 

V00;00(q) 



3aj_ 

x Erfc(|q|aj./V2) 
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can be rewritten as: 
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Note that V i ; oi = Viojoi = Vio ; io = Voi ; iq, V o ; ii = Vn ; oo 
and Vqo ; oi = Voo,io = Voijoo = lao;00 = 0. 



Appendix B: Numerical solution of the Hartree-Fock equation 



Our goal is to develop a numerical routine to find a self- 



consistent solution to Eq. ( 17 1 at any given temperature, dipo- 
lar interaction strength and transverse confinement strength. 



In its current form, Eq. ( 17 1 describes an infinite number of 
coupled non-linear integral equations. However, we note that 
at zero temperature, there is only a finite number of occupied 
subbands due to the sharp step-like behavior of Fermi occupa- 
tion function. At finite temperatures, we also expect to find a 
finite number of subbands with a significant population. Thus, 
a cut-off can be imposed on the number of subbands in prac- 
tice. In our implementation, we choose the cut-off such that 
the fractional density of the highest neglected subband is less 
than 1CP 3 . This ensures that the presence of higher subbands 
have a negligible effect, i.e. of the order of 10~ 3 Z? 2 hp, on 
the energies of the lower subbands (See Eq{T7|. The number 
of dimensions of the integral equations can also be reduced 
due to the isotropy of the inter-subband interactions, which is 
inherited by the self-energy matrices and their related quan- 
tities. The self-consistent Hartree-Fock self-energy equation 



xU(k') lp n F [£ p (k')}, (Bl) 

where W IJ ,x ] yi,(k, k') is the angle-averaged interaction, de- 
fined as: 

W^ n „(k, k')= J U V mA;7 „ (V<Z 2 + fc 2 -2gfcco S( />) . 

(B2) 

From a numerical perspective, it is favorable to deal with inte- 
grals with bounded integration domains. Although the mo- 
is unbounded, an up- 
per bound can be imposed on it in a controlled way. Since 
the subbands have (approximately) a quadratic energy disper- 
sion for large momenta, the Fermi occupation of states falls 
super-exponentially fast for large k' and the integral kernel of 
Eq. ( |B1| | becomes negligible for large k' . In our implemen- 
tation, we impose a large momentum cut-off A" cut , such that 
" F KpO')] < !0~ 6 for all k' > K cnt . Finally, we approxi- 
mate the resulting bounded integrations using quadrature for- 
mulas: 
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x Uik'^Uik'^y [U k 'j)] , i = l,..-,N p (B3) 

where fc, and w, denote the quadrature nodes and weights 
respectively and N p is the number of quadrature points. We 
used a 200-point Gauss-Lobatto quadrature in our implemen- 
tation eh. 

It is desirable from an experimental perspective to solve the 
Hartree-Fock equations for a given particle density n instead 
of the chemical potential fi. Therefore, the chemical potential 
must be found in a such a way that the self-energy equation 



(Eq. B3 1 and particle density equation, 



E 



d 2 k 
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are satisfied simultaneously. We found it more efficient to 



solve the system of non-linear equations (Eqs. |B3J and B4i 
using an interior reflective trust region method |62 1 instead of 
the usual solution by iterations. In all of the runs, the solution 
was unique and the convergence was rapid. 



Appendix C: Energy absorption rate for an exponentially 
switched-on a.c. modulation pulse 

In this appendix, we derive the formula given for the energy 
absorption rate for an exponentially switched-on a.c. modula- 



tion pulse in Sec. IV (Eq. 33 1. Our strategy is to evaluate the 
total absorbed energy for a given pulse shape, AE. We define 
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the energy absorption rate as AE/ At, where At is the effec- 
tive time interval during which the external field interacts with 
the system. For the exponentially switched on and off pulse 
of Eq. (32 1, At ~ 2rj~ 1 . Our starting point is Eq. (28 1, which 



we integrate to get AE: 

poo 

AE = te / dt/(t) 



■oo 

duj 
2^ 



dt'f(t')d t T(t-t') 



l/MI 2 w rM. 



(Cl) 



Here, f(ui) is the Fourier transform of Eq. (32 1. We write 

l/MI 2 as: 

| JM| 2 = fl H + f_ («) + 2/+(w)/_(w), (C2) 



where: 



J7 
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First, we focus on the contributions of to AE 1 . This 

function has two second order poles at frequencies =pfi + iij 
and ±57 — 17/. Since 7~(t) is a causal function, T(w) is analytic 
in the upper complex plane and the integral in the second line 



of Eq. (Cl I can be evaluated by closing the contour on the 
upper half-plane. The resulting contribution is easily found to 
be: 



AE 4 



i-2 2 



T(±£l + ir)) {±0. + irj)T (±& + irj) 



(2^)2 



2(±fi + ii])T(±n + irj) 
(2ir/) 3 
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where T'M = (d/dw)T(w). 



The contribution of the cross-term, 2/+(w)/_(a;), to AE 
can be evaluated in the same way and we get: 



AE c = l -^-[r(ii + iv)-r(-Q 



17]) 
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Using Eqs. (|29|l, (|30|l, (|38j and ([39j, it is straightforward to 
establish the following identities: 

T(lu) = T(-io), 

= T(u))*, 
T'{u) = -T'{-w), 

r'K) = r'(o;)*, (C6) 

using which we get T(— £7 + irj) = T{VL + irf)* and T'(— 57 + 
irj) = —T'{Q. + irj)* . Using these identities, we get the fol- 
lowing simple expression for E = r/AE/2 = rj (AE + + 
AE- + AE c )/2: 

e = | mt [T{n + irj) + (n + i^nn + irj)} 

- - 3m [{il + ir!)T{n + irj)] . (Cl) 

In the limit 77/Q <C 1, the term in the second line is dominant 
and we get the desired result (Eq.|33]l. 
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FIG. 11. The three relaxation processes contributing to 



Appendix D: Estimation of the relaxation rate of inter-subband 
p-h excitations 

In this appendix, we estimate the relaxation rate of inter- 
subband p-h excitations created by the lattice modulation 
pulse. As mentioned in the final remarks of Sec. IV A such 



relaxation processes are absent in the TDHF approximation. 
Therefore, predictions of TDHF approximation for the inter- 
subband transition rates are only valid if the lifetime of exci- 



tations are larger than the pulse duration (See Eq. 46 1. 

For weak interactions, we expect the relaxation rate of 
both types of inter-subband excitations, quasiparticle-like 
excitations and excitons, to be of the same order. Therefore, 
we study only the former case here. We are interested in 
relaxation processes that change the population of subbands 
and therefore, we neglect the intra-subband scatterings. 

Relaxation of a quasiparticle excitation can occur through 
several channels. The number of relaxation channels is 
higher for quasiparticles in higher subbands. As mentioned 
in Sec. |ll] the inter-subband interaction matrix elements 
conserve the net parity of the interacting particles and this 
condition constrains the number of relaxation channels. 
A simple combinatorial analysis shows that the number 
of relaxation channels scales like ~ N 3 /4, where N is 
the subband index of the decaying quasiparticle. One can, 
however, argue that the most important contributions still 
result from the processes involving only the first few lower 
subbands. The reason is twofold: (1) the inter-subband 
interaction between particles in well separated subbands is 
small, and (2) the lower subbands have larger Fermi momenta 
and consequently, creating a quasiparticle excitation in them 
upon de-exciting the original quasiparticle has a higher 
energy cost. Therefore, the quasiparticle relaxation rate of 
particles in the higher subbands are expected to be of the 
same order as those in the lower subbands. 
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FIG. 12. (a) the total momentum-averaged relaxation rate of a quasiparticle in the second excited subband, (r| ) as function of transverse 
confinement width aj_. The solid and dashed lines correspond to T = and T = 10 T* respectively, (b), (c) and (d) show the momentum- 
resolved contribution of different relaxation processes. The black and blue plots correspond to |po| = and |po| = pf,o respectively, where 
Pf,o = hy/4-Kn, is the Fermi momentum of the zeroth subband at zero temperature. The inset figures are graphical representations of the 
processes. 



As a concrete example, we study the relaxation of a quasi- 
particle excitation in the second excited subband, as a func- 
tion of transverse confinement width and temperature. Fig. 12 
shows the three relaxation processes that that contribute to the 
decay. Since the interactions are assumed to be weak (r^ <C 
1), Bom's approximation is applicable. We ignore the self- 
energy correction since such contributions are beyond Born's 
approximation. The initial and final wavefunctions of the sys- 
tem in all of the processes shown in Fig 
as \i) — cl ,c_ J'I'o) and If) = ct ~c, c , 
respectively, where l^o) is the Fermi gas-like ground state of 
the system. The rate of this generic process can be obtained 
using Fermi's golden rule: 



12 can be written 



r 7 A;^l/(PO,Po>Pl)Pl) — 
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We work in units in which h = 1 in this appendix. Since the 
particles are identical, T 1 \.^ l , = T^x-^ and over-counting 
must be carefully avoided. The total contribution of this pro- 
cess can be found by summing over p , pi and p[. After a 



lengthy but straightforward algebra, we obtain: 
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where Q = Q(cos0x + sin^y), q = <?o(cos tpx + sintpy) 
and q = y/(X + v - 7 - \i) 2mfej trap + |Q - p | 2 . The 
total relaxation rate is obtained by summing over all channels: 

r r 2 oh (Po) -r 02 . 00 (p ) + r 12;10 (p ) + r 02;11 ( Po ). (D3) 

We define the momentum-averaged relaxation rate as: 
d 2 p 



or- > = 



J (2tt) s 



[n (p)-n 2 ( P )] rf-(p) 
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[n (p) - n 2 (p) 



where the averaging is weighted according to the mo- 
mentum distribution of the p-h excitations generated 
by the lattice modulation pulse. For weak interactions, 
the number of excitations at momentum p is propor- 
tional to the equilibrium occupation of the states, i.e. 
oc n (p)[l - n 2 (p)] - n 2 (p)[l - n (p)] = n (p) - n 2 (p). 
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We evaluate the triple integral in Eq. ( D2 1 numerically us 
ing an adaptive Monte Carlo integration algorithm. Fig. [12} 



shows (r| ) as a function of transverse confinement width 



a ± for T = and T = 10 T*. Fig. [T2f)-d show the 
momentum-resolved rate of each of the relaxation processes 
as a function of a± for |po| = and |po| = pf.q, where 
Pf.q = h^/Aim, is the Fermi momentum of the zeroth sub- 
band at zero temperature. 

We find that the major contributions result from the first two 
processes, the former being dominant for stronger traps. Also, 
increasing the temperature and quasiparticle momenta natu- 
rally results in higher relaxation rates. According to Fig. [T2fr , 
the highest relaxation rate is (r| ) max ~ hr 2 d [Hn/2m] and 
10 r\ [hn/2m] for T = and T = 10 T* respectively. Us- 
ing these estimates and Eq. (46 1, choosing -q > 5 r 2 , [hn/2m] 



guarantees that the predictions of TDHF approximation for 
inter-subband transition rates are reliable for temperatures up 
toT = 10T*. Forr d = 0.05, we get 77 > 0.01 [hn/2m]. 



Appendix E: Numerical Evaluation of h a p(ijj) 



In this appendix, we discuss the numerical method used to 



solve Eq. (41 1. Our approach is similar to similar to that de- 
scribed in Appendix [B| i.e. we approximate the integral equa- 
tions using integral quadratures. However, given that the in- 
tegral equation to be solved is linear, the resulting system of 
equations is also linear. 

As a first step, we note that the dimension of the integral 
equations can be reduced due to the isotropy of the inter- 
subband interactions and Eq. ( |4"T) i reduces to: 

h a p(ki;uj) = T a/3 (fci) 

//(1& '<S( fc i' fc ')n(S(fc'; W )V(fc';^), 



2tt 
where: 



(El) 



2ji 



2tt 



2tt 



VuPwQs, k', Jk 2 + k' 2 - 2kk' cos fa) 



(E2) 



First, we study the subtler case of adiabatic switching (?; — » 0) 
for which u = O + iy — » £1 + i0 + . In this limit, it is useful to 
separate the regular and singular parts of fl^l [k' , u): 



where P.V. denotes the Cauchy's principal value integration. 
The delta function may also be written as: 

^(n + -&.(*)) =E , a (E4) 

where {qA is the set of solutions of fl + £p(q) — £ a (s) = 0. 
Inserting Eqs. (E3 1 and (E4 1 into Eq. ( |E1[ > and approximating 
the principal value integrals using quadratures, we get: 



A'" 



h aP (kf p ;u) = SaySpxT^k? ) - ^ E 

p,v j=l 

,(d+e)f U aP n u(kj )-Tl fJl (k j ) 



27T 



x V(fcf ;w) 



w + ^(*f)-W) 



EE 



WM )-t'M )l 



xwSJS(*^,^")v(9T;«) 



(E5) 



where {kf } and {wf } are quadrature nodes and weights 
for the principal value integrals involving h a p-,„, is the 
number of corresponding nodes and {q^} is the set of solu- 
tions of fi + i,(q) - ^(q) = 0. 

In the case of finite swtiching rate (77 > 0), u = £1 + irj 
and n^°](fc, Li) is regular on the real axis. Therefore, the last 
term in Eq. ( |E5| l arising from the poles on the real axis will be 
absent. 

In our implementation, we imposed a large momentum cut- 
off on the integrals based on the same criteria as described 
in Sec. [B] Quadratures for principal value integrations were 
generated using a combination of 16th-order Gauss-Lobatto 
quadrature in the proximity of the singular points and Simpson 
rule elsewhere ll6D . The resulting complex linear system of 
equations was solved for h a p (k" , u) and consequently, the 
function T(u>) was evaluated using the discretized version of 
Eq. |02b: 



7» 



E 

a/3 



olB 7 a.0 



2tt 



np{k£ P ) -n a {kf) 



j-h a (3, lX (kf p ;uj) 



1% ' 



u + ^{k)-Uk) 
m [n p {k) - n a (k)] 5(Q + £p(k) - £ a (k)). (E3) 



Again, the last term is absent for finite switching rates. 

Once h a p(to) and T(ui) are found, the enegry absorp- 
tion spectrum, the spectral weight of excitons and the 
inter-subband transition rates can be readily evaluated using 



20 



Eq. ( 33 1, (43 1 and ( |44) i, respectively. For finite pulse switching 
rates, the excitons will be broadened and Eq. p3j) yields the 



full spectrum, including the broadened excitons and Eq. ( 43 1 
is no longer needed. 
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